#include <emmintrin.h>
#include <iomanip>
#include <algorithm>
#include <vector>
#include <iterator>
#include <cmath>

#ifdef _OPENMP
#include <omp.h>
#endif

#ifdef __x86_64__
#define NOT_ME
#endif

/** \brief base for class RealTimer; class SysTimer; class UserTimer;
  \ingroup util
  */

#include "timer.h"

#define _M 1024

template<class T>
int add0(T * A, const T * B, const T * C,
	 const size_t &nb, const T&p)
{
	size_t i = 0 ;
	for ( ; i<nb ; ++i){
		A[i] = B[i] + C[i] ;
		if (A[i]>=p) A[i] -= p ;
	}
	return 0;
}


template<class T>
int add1(T * A, const T * B, const T * C,
	 const size_t &nb, const T &p)
{
	size_t i = nb ;
	for ( ; i-- ; ){
		A[i] = B[i] + C[i] ;
		if (A[i]>=p) A[i] -= p ;
	}
	return 0;
}

template<class T>
int add2(T * A, const T* B, const T* C,
	 const size_t &nb, const T &p)
{
	T tmp;
	size_t i = nb ;
	for ( ; i-- ; ++A,++B,++C){
		tmp = *B + *C ;
		if (tmp<p) ; else tmp -= p ;
		*A = tmp;
	}
	return 0;
}

template<class T>
int add10( std::vector<T> & A,
	   const std::vector<T> & B,
	   const std::vector<T> & C,
	   const size_t & nb, const T & p)
{
	for (size_t i = 0 ; i < nb ; ++i){
		A[i] = B[i] + C[i] ;
		if (A[i] < p) ; else A[i] -= p ;
	}

	return 0;
}

template<class T>
int add11( std::vector<T> & A,
	   const std::vector<T> & B,
	   const std::vector<T> & C,
	   const size_t & nb, const T & p)
{
	size_t i = nb ;
	for (;i--;){
		A[i] = B[i] + C[i] ;
		if (A[i] < p) ; else A[i] -= p ;
	}
	return 0;

}


#ifdef _OPENMP
template<class T>
int add0_par(T * A, const T * B, const T * C,
	     const size_t &nb, const T&p)
{
	size_t i;
#pragma omp parallel for
	for ( i = 0 ; i < nb ; ++i){
		A[i] = B[i] + C[i] ;
		if (A[i]<p) ; else A[i] -= p ;
	}
	return 0;
}

template<class T>
int add1_par(T * A, const T * B, const T * C,
	     const size_t &nb, const T &p)
{
	size_t k = 128 ;
	size_t M = nb/k ;
	size_t i,j;
#pragma omp parallel for
	for ( i = 0 ; i < k ; ++i)
	{
		for ( j = 0 ; j < M ; ++j) {
			A[i*M+j] = B[i*M+j] + C[i*M+j] ;
			if (A[i*M+j]>=p) A[i*M+j] -= p ;
		}
	}

	for ( j = M*k ; j < nb ; ++j)
	{
		A[j] = B[j] + C[j] ;
		if (A[j]>=p) A[j] -= p ;
	}


	return 0;
}

template<class T>
int add2_par(T * A, const T* B, const T* C,
	     const size_t nb, const T & p)
{
	size_t i ;
#pragma omp parallel for //schedule(static,512)
	for ( i = 0 ; i < nb ; ++i){
		T tmp = B[i] + C[i] ;
		if (tmp<p) ; else tmp -= p ;
		A[i] = tmp;
	}
	return 0;
}

template<class T>
int add3_par(T * A, const T* B, const T* C,
	     const size_t &nb, const T & p)
{
	size_t i ;
#pragma omp parallel for //schedule(static,512)
	for ( i = 0 ; i < nb ; ++i){
		T tmp = B[i] + C[i] ;
		if (tmp<p) ; else tmp -= p ;
		A[i] = tmp;
	}
	return 0;
}
#ifdef NOT_ME
int add_sse_par(float * A, const float * B, const float * C,
		const size_t &nb)
{
	size_t j = 0 ;
	__m128 av, bv, cv;

	for (j = 0; j < nb-3; j+=4,A+=4,B+=4,C+=4){
		av = _mm_load_ps(A) ;
		bv = _mm_load_ps(B) ;
		cv = _mm_load_ps(C) ;
		av = _mm_add_ps(bv, cv);
		_mm_store_ps(A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;

	return 0;
}

int add_sse_par(float * A, const float * B, const float * C,
		const size_t &nb, const float &p)
{
	size_t j = 0 ;
	__m128 av, bv, cv;

	for (j = 0; j < nb-3; j+=4,B+=4,C+=4){
		av = _mm_load_ps(A) ;
		bv = _mm_load_ps(B) ;
		cv = _mm_load_ps(C) ;
		av = _mm_add_ps(bv, cv);
		_mm_store_ps(A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}

	return 0 ;
}

int add_sse_par(double* A, const double* B, const double* C,
		const size_t &nb)
{
	size_t j = 0 ;
	__m128d av, bv, cv;
	for (j = 0; j < nb-1 ;j+=2,A+=2,B+=2,C+=2){
		av = _mm_load_pd(A) ;
		bv = _mm_load_pd(B) ;
		cv = _mm_load_pd(C) ;
		av = _mm_add_pd(bv, cv);
		_mm_store_pd(A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;
	//        j *= 2 ;
	//        for (; j < nb ; j++)
	//                *(A+j) = *(C+j) + *(B+j) ;

	return 0 ;
}

int add_sse_par(double* A, const double* B, const double* C,
		const size_t &nb, const double &p)
{
	size_t j = 0 ;
	__m128d av, bv, cv;
	for (j = 0; j < nb-1 ;j+=2,B+=2,C+=2){
		av = _mm_load_pd(A) ;
		bv = _mm_load_pd(B) ;
		cv = _mm_load_pd(C) ;
		av = _mm_add_pd(bv, cv);
		_mm_store_pd(A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}

	return 0 ;
}

#if 1
int add_sse_par(int * A, const int * B, const int * C,
		const size_t &nb, const double & p)
{
	size_t j = 0 ;
	__m128i av, bv, cv;
	for (j = 0; j < nb-3 ;j+=4,B+=4,C+=4){
		av = _mm_load_si128((__m128i*)A) ;
		bv = _mm_load_si128((__m128i*)B) ;
		cv = _mm_load_si128((__m128i*)C) ;
		av = _mm_add_epi64(bv, cv);
		_mm_store_si128((__m128i*)A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;

	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	return 0;
}
int add_sse_par(int * A, const int * B, const int * C,
		const size_t &nb)
{
	size_t j = 0 ;
	__m128i av, bv, cv;
	for (j = 0; j < nb-3 ;j+=4,A+=4,B+=4,C+=4){
		av = _mm_load_si128((__m128i*)A) ;
		bv = _mm_load_si128((__m128i*)B) ;
		cv = _mm_load_si128((__m128i*)C) ;
		av = _mm_add_epi64(bv, cv);
		_mm_store_si128((__m128i*)A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;
	return 0;
}
#endif

#endif
#endif

template<class T>
int add3(T * A, const T* B, const T* C,
	 const size_t nb, const T &p)
{
	size_t i = nb ;
	for ( ; i-- ; ++A,++B,++C){
		*A = *B + *C ;
		if (*A>=p)  *A -= p ;
	}
	return 0;
}

#if 0
template<class T>
int add3(T* A, const T* B, const T* C, const size_t nb, T p)
{
	T* Av = A ;
	const T * Bv = B;
	const T * Cv = C;
	for (size_t b = 0 ; b<nb/_M; ++b){
		for (size_t i =0 ; i < _M ; ++i,++Av,++Bv,++Cv){
			*Av = *Bv + *Cv ;
			if (*Av>=p) *Av -= p;
		}
	}
	for ( ; Av < A+nb ; ++Av,++Bv,++Cv) {
		*Av = *Bv + *Cv ;
		if (*Av>=p) *Av -= p;
	}

	return 0;
}
#endif
#ifdef NOT_ME
int add_sse(float * A, const float * B, const float * C,
	    const size_t &nb)
{
	size_t j = 0 ;
	__m128 av, bv, cv;

	for (j = 0; j < nb-3; j+=4,A+=4,B+=4,C+=4){
		av = _mm_load_ps(A) ;
		bv = _mm_load_ps(B) ;
		cv = _mm_load_ps(C) ;
		av = _mm_add_ps(bv, cv);
		_mm_store_ps(A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;

	return 0;
}

int add_sse(float * A, const float * B, const float * C,
	    const size_t &nb, const float &p)
{
	size_t j = 0 ;
	__m128 av, bv, cv;

	for (j = 0; j < nb-3; j+=4,B+=4,C+=4){
		av = _mm_load_ps(A) ;
		bv = _mm_load_ps(B) ;
		cv = _mm_load_ps(C) ;
		av = _mm_add_ps(bv, cv);
		_mm_store_ps(A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}

	return 0 ;
}

int add_sse(double* A, const double* B, const double* C,
	    const size_t &nb)
{
	size_t j = 0 ;
	__m128d av, bv, cv;
	for (j = 0; j < nb-1 ;j+=2,A+=2,B+=2,C+=2){
		av = _mm_load_pd(A) ;
		bv = _mm_load_pd(B) ;
		cv = _mm_load_pd(C) ;
		av = _mm_add_pd(bv, cv);
		_mm_store_pd(A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;
	//        j *= 2 ;
	//        for (; j < nb ; j++)
	//                *(A+j) = *(C+j) + *(B+j) ;

	return 0 ;
}

int add_sse(double* A, const double* B, const double* C,
	    const size_t &nb, const double &p)
{
	size_t j = 0 ;
	__m128d av, bv, cv;
	for (j = 0; j < nb-1 ;j+=2,B+=2,C+=2){
		av = _mm_load_pd(A) ;
		bv = _mm_load_pd(B) ;
		cv = _mm_load_pd(C) ;
		av = _mm_add_pd(bv, cv);
		_mm_store_pd(A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}

	return 0 ;
}

#if 1
int add_sse(int * A, const int * B, const int * C,
	    const size_t &nb, const double &p)
{
	size_t j = 0 ;
	__m128i av, bv, cv;
	for (j = 0; j < nb-3 ;j+=4,B+=4,C+=4){
		av = _mm_load_si128((__m128i*)A) ;
		bv = _mm_load_si128((__m128i*)B) ;
		cv = _mm_load_si128((__m128i*)C) ;
		av = _mm_add_epi64(bv, cv);
		_mm_store_si128((__m128i*)A,av);
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;
		if (*A<p) ++A ; else *(A++) -= p ;

	}
	for (; j < nb ; ++j,++B,++C){
		*A = *C + *B ;
		if (*A<p) ++A ; else *(A++) -= p ;
	}
	return 0;
}
int add_sse(int * A, const int * B, const int * C,
	    const size_t &nb)
{
	size_t j = 0 ;
	__m128i av, bv, cv;
	for (j = 0; j < nb-3 ;j+=4,A+=4,B+=4,C+=4){
		av = _mm_load_si128((__m128i*)A) ;
		bv = _mm_load_si128((__m128i*)B) ;
		cv = _mm_load_si128((__m128i*)C) ;
		av = _mm_add_epi64(bv, cv);
		_mm_store_si128((__m128i*)A,av);
	}
	for (; j < nb ; ++j,++A,++B,++C)
		*A = *C + *B ;
	return 0;
}
#endif
#endif
#ifndef _OPENMP
template<class T>
void reduce(T*A,
	    const size_t &m,const  size_t &n,const  size_t &lda,const  T &p)
{
	if (n==lda){
		T* Aend = A+m*n;
		for ( ;  A<Aend ; ++A) if (*A>=p) *A -= p;
		// :) //
	}else{
		size_t i ;
		for (i=0 ; i<m;++i)
			for (size_t j = 0 ; j < n ; ++j)
				if (A[i*lda+j]>=p) A[i*lda+j] -= p;
	}
}
#else
template<class T>
void reduce(T*A,
	    const size_t &m,const  size_t &n,const  size_t &lda,const  T &p)
{
	size_t i = 0 ;
#pragma omp parallel for
	for (i=0 ; i<m;++i){
		T* Ai = A+i*lda;
		for (size_t j = 0 ; j < n ; ++j)
			if (Ai[j]<p) ; else Ai[j] -= p;
	}
}
#endif

template<class T>
void mul0(T*A,const T*B,const T*C,
	  const size_t &m, const size_t &n, const size_t &lda, const  T &p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add0_par(A,B,C,mn,p);
#else
		add0(A,B,C,mn,p);
#endif
	}else{
		size_t i =0 ;
#ifdef _OPENMP
#pragma omp parallel for
		for ( i=0; i < m ; ++i)
			add0(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		for ( i=0; i < m ; ++i, A+=lda,C+=lda,B+=lda)
			add0(A,B,C,n,p);
#endif
	}
	return;
}

template<class T>
void mul1(T*A,const T*B,const T*C,
	  const size_t &m, const size_t &n, const size_t &lda, const T &p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add1_par(A,B,C,mn,p);
#else
		add1(A,B,C,mn,p);
#endif
	}else{
		size_t i =0 ;
#ifdef _OPENMP
#pragma omp parallel for
		for (i=0 ; i < m ; ++i)
			add1(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		for (i=0 ; i < m ; ++i, A+=lda,C+=lda,B+=lda)
			add1(A,B,C,n,p);
#endif
	}
	return;
}

template<class T>
void mul2(T*A,const T*B, const T*C,
	  const size_t & m, const size_t &n, const size_t &lda, const T &p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add2_par(A,B,C,mn,p);
#else
		add2(A,B,C,mn,p);
#endif
	}else{

#ifdef _OPENMP
		size_t i;
#pragma omp parallel for
		for (i = 0 ; i < m ; ++i)
			add2(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		__restrict__ T* Aloc = A ;
		const __restrict__ T* Bloc = B ;
		const __restrict__ T* Cloc = C ;

		for (size_t i = 0 ; i < m ; ++i, Cloc+=lda,Aloc+=lda,Bloc+=lda)
			add2(Aloc,Bloc,Cloc,n,p);
#endif
	}
	return;
}

template<class T>
void mul3(T*A,const T*B,const T*C,
	  const size_t &m, const size_t &n, const size_t &lda, const T &p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add3_par(A,B,C,mn,p);
#else
		add3(A,B,C,mn,p);
#endif
	}else{

#ifdef _OPENMP
		size_t i;
#pragma omp parallel for
		for (i = 0 ; i < m ; ++i)
			add3(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		T* Aloc = A ;
		const T* Bloc = B ;
		const T* Cloc = C ;

		for (size_t i = 0 ; i < m ; ++i, Cloc+=lda,Aloc+=lda,Bloc+=lda)
			add3(Aloc,Bloc,Cloc,n,p);
#endif
	}
	return;
}

#ifdef NOT_ME
template<class T>
void mul4(T*A,const T*B,const T*C,
	  const size_t &m, const size_t &n, const size_t &lda, const T &p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add_sse_par(A,B,C,mn);
#else
		add_sse(A,B,C,mn);
#endif
		reduce(A,m,n,lda,p);
	}else{

#ifdef _OPENMP
		size_t  i;
#pragma omp parallel for
		for (i = 0 ; i < m ; ++i)
			add_sse(A+i*lda,B+i*lda,C+i*lda,n);

#else
		T* Aloc = A ;
		const T* Bloc = B ;
		const T* Cloc = C ;

		for (size_t i = 0 ; i < m ; ++i, Aloc+=lda,Bloc+=lda,Cloc+=lda)
			add_sse(Aloc,Bloc,Cloc,n);
#endif
		reduce(A,m,n,lda,p);
		return;
	}
}

template<class T>
void mul5(T*A,const T*B,const T*C,
	  const size_t& m, const size_t& n, const size_t& lda, const T & p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add_sse_par(A,B,C,mn,p);
#else
		add_sse(A,B,C,mn,p);
#endif
	}else{


#ifdef _OPENMP
		size_t i ;
#pragma omp parallel for
		for ( i = 0 ; i < m ; ++i)
			add_sse(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		T* Aloc = A ;
		const T* Bloc = B ;
		const T* Cloc = C ;

		for (size_t i = 0 ; i < m ; ++i, Aloc+=lda,Bloc+=lda,Cloc+=lda)
			add_sse(Aloc,Bloc,Cloc,n,p);
#endif
	}
	return;
}
#endif

template<class T>
void mul6(T*A,const T*B,const T*C,
	  const size_t& m, const size_t& n, const size_t& lda,const  T& p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add3_par(A,B,C,mn,p);
#else
		add3(A,B,C,mn,p);
#endif
	}else{


#ifdef _OPENMP
		size_t i =0 ;
#pragma omp parallel for
		for (i=0 ; i<m  ; ++i)
			add3(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		T* Aloc = A+lda*(m-1) ;
		const T* Bloc = B+lda*(m-1) ;
		const T* Cloc = C+lda*(m-1) ;
		size_t i=m ;

		for ( ; i--  ; Cloc-=lda,Aloc-=lda,Bloc-=lda)
			add3(Aloc,Bloc,Cloc,n,p);
#endif
	}
	return;
}

template<class T>
void mul7(T*A,const T*B,const T*C,
	  const size_t& m, const size_t& n, const size_t& lda,const  T& p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add3_par(A,B,C,mn,p);
#else
		add3(A,B,C,mn,p);
#endif
	}else{


#ifdef _OPENMP
		size_t i;
#pragma omp parallel for
		for (i =0  ; i<m  ; ++i)
			add3(A+i*lda,B+i*lda,C+i*lda,n,p);
#else

		T* Aloc = A ;
		const T* Bloc = B ;
		const T* Cloc = C ;
		size_t i=m ;

		for ( ; i--  ; )
			add3(Aloc+i*lda,Bloc+i*lda,Cloc+i*lda,n,p);
#endif
	}
	return;
}

template<class T>
void mul8(T*A,const T*B,const T*C,
	  const size_t& m, const size_t& n, const size_t& lda,const  T& p)
{
	if (n==lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add2_par(A,B,C,mn,p) ;
#else
		add2(A,B,C,mn,p) ;
#endif
	}else {

#ifdef _OPENMP
		size_t i = 0 ;
#pragma omp parallel for
		for ( i=0; i<m  ;++i)
			add2(A+i*lda,B+i*lda,C+i*lda,n,p);
#else
		size_t i=m ;

		for ( ; i--  ; C+=lda,A+=lda,B+=lda)
			add2(A,B,C,n,p);
#endif
	}
	return;
}

#ifdef NOT_ME
template<class T>
void mul9(T*A,const T*B,const T*C,
	  const size_t& m, const size_t& n, const size_t& lda,const  T& p)
{
	if (n == lda) {
		size_t mn = m*n ;
#ifdef _OPENMP
		add_sse_par(A,B,C,mn,p);
#else
		add_sse(A,B,C,mn,p);
#endif
	}else{


#ifdef _OPENMP
		size_t i;
#pragma omp parallel for
		for (i = 0 ; i<m ; ++i)
			add_sse(A+i*lda,B+i*lda,C+i*lda,n,p);

#else
		size_t i=m ;
		for (; i-- ; A+=lda,B+=lda,C+=lda)
			add_sse(A,B,C,n,p);
#endif
	}
	return;
}
#endif


template<class T>
void mul10(std::vector<std::vector<T> > & A,
	   const std::vector<std::vector<T> > & B,
	   const std::vector<std::vector<T> > & C,
	   const size_t &m, const size_t &n, const T &p)
{
	size_t i ;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (i = 0 ; i < m ; ++i)
		add10(A[i],B[i],C[i],n,p);
}

template<class T>
void mul11(std::vector<std::vector<T> > & A,
	   const std::vector<std::vector<T> > & B,
	   const std::vector<std::vector<T> > & C,
	   const size_t &m, const size_t &n, const T &p)
{
	size_t i ;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (i = 0 ; i < m ; ++i)
		add11(A[i],B[i],C[i],n,p);
}

template<class T>
bool correct(const T*A, const T*B, const T*C,
	     const size_t &m, const size_t &n, const size_t & lda, const T&p)
{
	bool ret = true ;
	for (size_t i = 0 ; i < m && ret ; ++i)
		for (size_t j = 0 ; j < n && ret ; ++j)
			if (A[i*lda+j] != (fmod((double)B[i*lda+j]+C[i*lda+j],(double)p)))
				ret = false ;
	return ret ;
}
template<class T>
bool correct(const std::vector<std::vector<T> > &A,
	     const std::vector<std::vector<T> > &B,
	     const std::vector<std::vector<T> > &C,
	     const size_t &m, const size_t &n, const T&p)
{
	bool ret = true ;
	for (size_t i = 0 ; i < m && ret ; ++i)
		for (size_t j = 0 ; j < n && ret ; ++j)
			if (A[i][j] != (fmod((double)B[i][j]+C[i][j],(double)p)))
				ret = false ;
	return ret ;
}



template<class T>
void launch(const size_t & m, const size_t &n, const size_t &lda, const size_t &M, const T &p)
{
	std::cout << "prêts ?" << std::endl;
	//size_t lda = 16*(size_t)std::ceil((double)n/(double)16); // ALIGNED !!
	size_t mn = m*lda ;
	//        std::cout << mn << std::endl;
	T *a = new T[mn];
	T *b = new T[mn];
	T *c = new T[mn];
	const char * comments[10] ;

	const std::vector<T> Z(lda,0);
	std::vector<std::vector<T> > A(m,Z), B(m,Z),C(m,Z);
	for (size_t i = 0 ; i < mn ; ++i)
		b[i] = c[i] = fmod(i,p) ;
	for (size_t i = 0 ; i < m ; ++i)
		for (size_t j = 0 ; j < n ; ++j)
			B[i][j] = C[i][j] = fmod(i*lda+j,p) ;

#ifdef NOT_ME
	size_t N = 10 ;
#else
	size_t N = 7 ;
#endif
	std::cout << "partez !" << std::endl;

	if (M != 1)
	{
		std::vector<double> TT(N,0) ;
#ifdef _OPENMP
		OMPTimer Tim;
#else
		Timer Tim ;
#endif
		size_t k = 0 ;
		/*  Timings */
		comments[k]="i++  A[i] <" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul0(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!(correct<T>(a,b,c,m,n,lda,p))) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

		comments[k]="i--  A[i] <" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul1(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

		comments[k]="i--  *A    <" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul2(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

		comments[k]="i--  *A    >=" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul3(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

#ifdef NOT_ME
#if 0
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul4(a,b,c,m,n,lda,p);
		Tim.stop() ; TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
#endif

		comments[k]="sse" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul5(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());
#endif

		comments[k]="i--  *A    >=  j--  A-=" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul6(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

		comments[k]="i--  *A    >=  j--  A+j" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul7(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

		comments[k]="i--  *A    <   j--  A+=" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul8(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

#ifdef NOT_ME
		comments[k]="sse i--" ;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul9(a,b,c,m,n,lda,p);
		Tim.stop() ;
		if (!correct<T>(a,b,c,m,n,lda,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());
#endif

		comments[k]="vector i++" ;
		for (size_t i = 0 ; i < m ; ++i) for (size_t j = 0 ; j < n; ++j) A[i][j] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul10(A,B,C,m,n,p);
		Tim.stop() ;
		if (!correct<T>(A,B,C,m,n,p)) std::cout << "* erreur " << k << std::endl;
		TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());

#if 0
		comments[k]="vector i--" ;
		for (size_t i = 0 ; i < m ; ++i) for (size_t j = 0 ; j < n; ++j) A[i][j] =0 ;
		Tim.clear() ; Tim.start() ;
		for (size_t i = 0 ; i < M; ++i)
			mul11(A,B,C,m,n,p);
		Tim.stop() ; TT[k++] = (double)M*n*m/(1.e6*Tim.usertime());
		if (!correct<T>(A,B,C,m,n,p)) std::cout << "* erreur " << k << std::endl;
#endif



		typedef 	std::vector<double>::iterator iter ;
		iter it = std::max_element(TT.begin(),TT.end());
		size_t i = 0 ;
		for (iter jt = TT.begin(); jt != TT.end() ; ++jt,++i ) {
			std::cout << "mul" << i << "\t\033[" ;
			if (std::abs(*jt-*it)/(*it)<0.02)
				std::cout << "1;32m";
			else if (std::abs(*jt-*it)/(*it)<0.05)
				std::cout << "0;32m";
			else if (std::abs(*jt-*it)/(*it)<0.1)
				std::cout << "0;33m";
			else if (std::abs(*jt-*it)/(*it)<0.15)
				std::cout << "0;35m";
			else if (std::abs(*jt-*it)/(*it)<0.2)
				std::cout << "0;31m";
			else
				std::cout << "1;31m";
			std::cout <<*jt << "\033[0m\tMflops\t|" << comments[i] << std::endl;
		}


	}
	else
	{
		std::cout << "mul0" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul0(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;


		std::cout << "mul1" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul1(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

		std::cout << "mul2" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul2(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

		std::cout << "mul3" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul3(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

#ifdef NOT_ME
#if 0
		std::cout << "mul4" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul4(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;
#endif

		std::cout << "mul5" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul5(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;
#endif

		std::cout << "mul6" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul6(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

		std::cout << "mul7" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul7(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

		std::cout << "mul8" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul8(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;

#ifdef NOT_ME
		std::cout << "mul9" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul9(a,b,c,m,n,lda,p);
		for (size_t i = 0 ; i < mn ; ++i) std::cout << a[i] << ' ' ; std::cout << std::endl;
#endif
		std::cout << "mul10" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul10(A,B,C,m,n,p);
		for (size_t i = 0 ; i < m ; ++i) for (size_t j = 0 ; j < lda ; ++j) std::cout << A[i][j] << ' ' ; std::cout << std::endl;

		std::cout << "mul11" << std::endl;
		for (size_t i = 0 ; i < mn ; ++i) a[i] =0 ;
		mul11(A,B,C,m,n,p);
		for (size_t i = 0 ; i < m ; ++i) for (size_t j = 0 ; j < lda ; ++j) std::cout << A[i][j] << ' ' ; std::cout << std::endl;




	}

	delete [] a;
	delete [] b;
	delete [] c;
}

int main(int ac, char ** av)
{
	int M = ac>1?atoi(av[1]):10;
	int m = ac>2?atoi(av[2]):1024;
	int n = ac>3?atoi(av[3]):1024;
	int p = ac>4?atoi(av[4]):113;
	size_t lda = 16*(size_t)std::ceil((double)n/(double)16); // ALIGNED !!
	n = lda;

	std::cout << "répet " << M << "x matrices " << m << 'x' << n << " (" << lda <<')' << std::endl;
	std::cout << " * f" << std::endl;
	launch<float>(m,n,lda,M,p) ;
	std::cout << " * d" << std::endl;
	launch<double>(m,n,lda,M,p);
	std::cout << " * i" << std::endl;
	launch<int>(m,n,lda,M,p);

	if (M!=1) lda += 1024 ;
	else lda +=16 ;
	std::cout << "répet " << M << "x matrices " << m << 'x' << n << " (" << lda <<')' << std::endl;
	std::cout << " * f" << std::endl;
	launch<float>(m,n,lda,M,p) ;
	std::cout << " * d" << std::endl;
	launch<double>(m,n,lda,M,p);
	std::cout << " * i" << std::endl;
	launch<int>(m,n,lda,M,p);

	return 0;
}
/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s:syntax=cpp.doxygen:foldmethod=syntax
